Modern digital cameras use an electronic sensor for acquiring the spatial variations of light intensity. The two main fabrication technologies of sensors are the so called CCD (Charge Coupled Device) and the CMOS (Complementary Metal Oxide Semiconductor).
For example, a CCD is provided by an array of photosensitive diodes that converts photons into electrons (electric charge). Each diode, called a photodiode produces an electric charge proportional to the quantity of light reaching it.
The single sensors of a digital camera are normally arranged as a two-dimensional array of millions of cells (photodiodes), which convert an infinitesimal portion of the image into electric charge. The stored generated charge is then transformed into a pixel of digital image by an analog-to-digital converter.
A monochromatic sensor is made sensitive to a single color component by a color filter. A filter may act on primary colors or on complementary colors. The primary color filters are commonly arranged according to a scheme known in the art by the name of Bayer pattern in which half of the pixels are made sensitive to the green color, a fourth of the pixels to the red color and the last fourth to the blue color. As depicted in FIG. 1, green and blue pixels alternate along a row, while green and red pixels alternate on the successive row and this scheme repeats itself for the whole image.
The main reason the number of acquired green pixels are doubled compared to the number of acquired red and blue pixels is the greater sensitivity of human vision of the green chromatic component. Therefore, the value of each pixel of a CCD sensor carries information on a single primary color.
To reconstruct a good quality color image it is therefore desired to carry out an operation known as color interpolation for generating a triplet of values RGB for each pixel of the image to be produced at an output. Therefore, for example, for estimating the red component of a pixel for which only the green component has been acquired by the filtered light sensitive cell, it is necessary to take into consideration the adjacent pixels carrying information on the red component.
Bicubic and bilinear color interpolation are simple and efficient algorithms, but they may introduce large errors in the contour zones of the images. To address this problem several authors have proposed optimization algorithms that are sensitive to the input data. For example, Adams introduced an edge oriented method, and Kimmel, Adams, and Pei-Tam suggested other color correlation methods.
Triangulation is another well-known technique of linear interpolation that has been proposed for performing zooming on a digital image. A brief review of triangulation-based linear interpolation approaches is presented herein below.
Triangulation: Basic Considerations
Basically, the approach is to partition an image into triangles, the vertices of which are the known pixels, according to certain criteria. The color rendering or reconstruction phase follows during which each pixel value of the image to be produced is calculated through a linear interpolation of the vertices of the pertinent triangle.
Given a set of distinct points V={(xi, yi)}, what is of interest is a convex hull triangulation of V, that is a set T={Ti}i=1 . . . q of non degenerate triangles that satisfies the following conditions:    every triangle vertex is an element of V, and every element of V is a triangle vertex;    every edge of a triangle in T contains exactly two points from V;    the intersection of any two different triangles in T is either empty, or is a shared edge or a shared vertex; and    the union of all triangles in T is the convex hull of V.
The criteria for choosing one among all the possible triangulations is to minimize a certain cost function C(T) of the entire triangulation T. Dyn-Levin-Rippa [5] suggested, among others, the following overall function cost:C(T)=Σ|cost(e)|  (1)for each edge e of T with a chosen cost(e). Once a possible triangulation is chosen, the values of the pixels will be interpolated by the unique 3D linear polynomial passing for the three vertices of the triangles. The 3D liner polynomial will be:Pi(x,y)=aix+biy+ci  (2)for each triangle Ti of T.Delaunay Triangulation
According to the Delaunay triangulation each triangle is such that the circle passing for his three vertices does not contain any other vertex of V, where V={vi} set of distinct vertices. The construction of a Delaunay triangulation is defined through its Voronoi diagram.
The Voronoi diagram for V is the partition Vor(V) of the plane into the Voronoi polygons associated to V. The Voronoi polygon associated to vi is the locus of points in the plane that are closer to vi than to any other member of V.
The Delaunay triangulation can be defined as follow:    partitioning of the set V in two subsets V1 and V2 of the same dimension;    recursively constructing Vor(V1) and Vor(V2); and    merging of Vor(V1) and Vor(V2) to obtain Vor(V).FIGS. 2 and 3 show the duality of the two problems.
If the vertices of V are equally spaced (as it would be in the case of color interpolation), the resulting triangles will be equilateral triangles. Dyn-Levin-Rippa [9] demonstrated that triangulations involving long, narrow triangles give better results in linear interpolation.
Data Dependent Triangulation
Several optimization algorithms for performing what may be defined as Data Dependent Triangulation have been developed. For example, Lawson's algorithm [10] converges to the globally optimal solution if the cost function is based on the Delaunay criterion. Schumaker's simulated annealing method [11], gives good results, but is relatively slow. The “edge swap with look-ahead” method [8], is relatively faster than other algorithms.
The Yu-Morse-Sederberg's “edge swap with look-ahead” method includes the following steps:                a) the phase of triangle construction (step 0) is basically a Delaunay Triangulation, conducted in a very simple way: all the vertices are joined to form quadrilaterals (for vertices equally spaced they would be squares) and the diagonal (edge) having the lower cost according to the particular cost function chosen in the quadrilateral will determine the initial triangulation, as depicted in FIG. 4. FIG. 5 illustrates a single diagonal swapping inside a considered quadrilateral; and        b) an iterative phase (step 1 . . . n) during which each quadrilateral and its four adjacent triangles are considered, trying all the possible single swapping of the diagonals to minimize the local cost (that is the cost of all the 13 edges involved). The routine ends when all quadrilaterals have been processed. FIG. 6 illustrates, at bottom, the possible 4 triangulations originated from a starting triangulation (depicted at the top of the figure) by a single edge swap.The process is repeated n times, using the results of the previous iteration step.Cost Functions        
A number of cost functions for characterizing a satisfactory triangulation interpolation have been proposed.
In particular, Dyn-Levin-Rippa [5] identify four cost functions called NC1 (Nearly C1), as particularly efficient. The cost for an edge e=[(x1,y1), (x2,y2)] in common with two triangles T1 and T2, the interpolating planes of which are respectively P1(x,y)=a1x+b1y+c1 and P2(x,y)=a2x+b2y+c2, may be calculated according to a selected cost function:    a) ABN (Angle Between Normals) cost function disclosed in Choi, et al. [12]cost(e)=cos(α)  (3)            where α is the angle between the normals to the interpolating planes P1 and P2.            b) JND (Jump in Normal Derivatives) cost functioncost(e)=|nx(a1−a2)+ny(b1−b2)|  (4)            where (nx, ny) is a unit vector orthogonal to the projection of the edge e.            c) DLP (Deviation from Linear Polynomials) cost functioncost(e)=∥h∥  (5)            where        
      h    =          [                                                                                                                      P                    1                                    ⁡                                      (                                                                  x                        2                                            ,                                              y                        2                                                              )                                                  -                                  F                  2                                                                                                                                                                                          P                    2                                    ⁡                                      (                                                                  x                        1                                            ,                                              y                        1                                                              )                                                  -                                  F                  1                                                                                      ]        ;                 Fi=value (xi, yi)            d) DP (Distance from Planes) cost functioncost(e)=∥g∥  (6)            where        
      g    =          [                                                  dist              ⁡                              (                                                      P                    1                                    ,                                      w                    2                                                  )                                                                                        dist              ⁡                              (                                                      P                    2                                    ,                                      w                    1                                                  )                                                        ]        ;                 wi=(xi, yi, Fi);        
      dist    ⁡          (                        P          j                ,                  w          t                    )        =                                                P            j                    ⁡                      (                                          x                t                            ,                              y                t                                      )                          -                  F          t                                  sqrt      ⁡              (                              a            j            2                    +                      b            j            2                    +          1                )            
Another cost function has been proposed by Yu-Morse-Sederberg [8]:cost(e)=∥□P1∥*∥□P2∥−□P1*□P2  (7)where∥□Pi∥=sqrt(ai2+bi2)Pi(x,y)=aix+biy+ci□Pi=(ai, bi) gradient of Pi 